Quasi-Newton method

In optimization, quasi-Newton methods (also known as variable metric methods) are algorithms for finding local maxima and minima of functions. Quasi-Newton methods are based on Newton's method to find the stationary point of a function, where the gradient is 0. Newton's method assumes that the function can be locally approximated as a quadratic in the region around the optimum, and use the first and second derivatives to find the stationary point.

In Quasi-Newton methods the Hessian matrix of second derivatives of the function to be minimized does not need to be computed. The Hessian is updated by analyzing successive gradient vectors instead. Quasi-Newton methods are a generalization of the secant method to find the root of the first derivative for multidimensional problems. In multi-dimensions the secant equation is under-determined, and quasi-Newton methods differ in how they constrain the solution, typically by adding a simple low-rank update to the current estimate of the Hessian.

The first quasi-Newton algorithm was proposed by W.C. Davidon, a physicist working at Argonne National Laboratory. He developed the first quasi-Newton algorithm in 1959: the DFP updating formula, which was later popularized by Fletcher and Powell in 1963, but is rarely used today. The most common quasi-Newton algorithms are currently the SR1 formula (for symmetric rank one), the BHHH method, the widespread BFGS method (suggested independently by Broyden, Fletcher, Goldfarb, and Shanno, in 1970), and its low-memory extension, L-BFGS. The Broyden's class is a linear combination of the DFP and BFGS methods.

The SR1 formula does not guarantee the update matrix to maintain positive-definiteness and can be used for indefinite problems. The Broyden's method does not require the update matrix to be symmetric and it is used to find the root of a general system of equations (rather than the gradient) by updating the Jacobian (rather than the Hessian).

1 Description of the method
2 Example algorithm using matlab
3 See also
4 References

Description of the method

As in Newton's method, one uses a second order approximation to find the minimum of a function $f(x)$ . The Taylor series of $f(x)$ around an iterate is:

$f(x_k%2B\Delta x) \approx f(x_k)%2B\nabla f(x_k)^T \Delta x%2B\frac{1}{2} \Delta x^T {B} \, \Delta x,$

where ( $\nabla f$ ) is the gradient and $B$ an approximation to the Hessian matrix. The gradient of this approximation (with respect to $\Delta x$ ) is

$\nabla f(x_k%2B\Delta x) \approx \nabla f(x_k)%2BB \, \Delta x$

and setting this gradient to zero provides the Newton step:

$\Delta x=-B^{-1}\nabla f(x_k), \,$

The Hessian approximation $B$ is chosen to satisfy

$\nabla f(x_k%2B\Delta x)=\nabla f(x_k)%2BB \, \Delta x,$

which is called the secant equation(the Taylor series of the gradient itself). In more than one dimension $B$ is under determined. In one dimension, solving for $B$ and applying the Newton's step with the updated value is equivalent to the secant method. Various methods are used to find the solution to the secant equation that is symmetric ( $B^T=B$ ) and closest to the current approximate value $B_k$ according to some metric $\min_B \|B-B_k\|$ . An approximate initial value of $B_0=I$ is often sufficient to achieve rapid convergence. The unknown $x_k$ is updated applying the Newton's step calculated using the current approximate Hessian matrix $B_{k}$

$\Delta x_k=- \alpha_k B_k^{-1}\nabla f(x_k)$ , with $\alpha$ chosen to satisfy the Wolfe conditions;
$x_{k%2B1}=x_{k}%2B\Delta x_k$ ;
The gradient computed at the new point $\nabla f(x_{k%2B1})$ , and

$y_k=\nabla f(x_{k%2B1})-\nabla f(x_k),$

is used to update the approximate Hessian $\displaystyle B_{k%2B1}$ , or directly its inverse $\displaystyle H_{k%2B1}=B_{k%2B1}^{-1}$ using the Sherman-Morrison formula.

A key property of the BFGS and DFP updates is that if $B_k$ is positive definite and $\alpha_k$ is chosen to satisfy the Wolfe conditions then $\displaystyle B_{k%2B1}$ is also positive definite.

The most popular update formulas are:

Method	$\displaystyle B_{k%2B1}=$	$H_{k%2B1}=B_{k%2B1}^{-1}=$
DFP	$\left (I-\frac {y_k \, \Delta x_k^T} {y_k^T \, \Delta x_k} \right ) B_k \left (I-\frac {\Delta x_k y_k^T} {y_k^T \, \Delta x_k} \right )%2B\frac{y_k y_k^T} {y_k^T \, \Delta x_k}$	$H_k %2B \frac {\Delta x_k \Delta x_k^T}{y_k^{T} \, \Delta x_k} - \frac {H_k y_k y_k^T H_k^T} {y_k^T H_k y_k}$
BFGS	$B_k %2B \frac {y_k y_k^T}{y_k^{T} \Delta x_k} - \frac {B_k \Delta x_k (B_k \Delta x_k)^T} {\Delta x_k^{T} B_k \, \Delta x_k}$	$\left (I-\frac {y_k \Delta x_k^T} {y_k^T \Delta x_k} \right )^T H_k \left (I-\frac { y_k \Delta x_k^T} {y_k^T \Delta x_k} \right )%2B\frac {\Delta x_k \Delta x_k^T} {y_k^T \, \Delta x_k}$
Broyden	$B_k%2B\frac {y_k-B_k \Delta x_k}{\Delta x_k^T \, \Delta x_k} \, \Delta x_k^T$	$H_{k}%2B\frac {(\Delta x_k-H_k y_k) y_k^T H_k}{y_k^T H_k \, \Delta x_k}$
Broyden family	$(1-\varphi_k) B_{k%2B1}^{BFGS}%2B \varphi_k B_{k%2B1}^{DFP}, \qquad \varphi\in[0,1]$
SR1	$B_{k}%2B\frac {(y_k-B_k \, \Delta x_k) (y_k-B_k \, \Delta x_k)^T}{(y_k-B_k \, \Delta x_k)^T \, \Delta x_k}$	$H_{k}%2B\frac {(\Delta x_k-H_k y_k) (\Delta x_k-H_k y_k)^T}{(\Delta x_k-H_k y_k)^T y_k}$

Example algorithm using matlab

%***********************************************************************%
% Usage: [x,Iter,FunEval,EF] = Quasi_Newton (fun,x0,MaxIter,epsg,epsx)
%         fun: name of the multidimensional scalar objective function
%              (string). This function takes a vector argument of length
%              n and returns a scalar.
%          x0: starting point (row vector of length n).
%     MaxIter: maximum number of iterations to find a solution.
%        epsg: maximum acceptable Euclidean norm of the gradient of the
%              objective function at the solution found.
%        epsx: minimum relative change in the optimization variables x.
%           x: solution found (row vector of length n).
%        Iter: number of iterations needed to find the solution.
%     FunEval: number of function evaluations needed.
%          EF: exit flag,
%              EF=1: successful optimization (gradient is small enough).
%              EF=2: algorithm converged (relative change in x is small 
%                    enough).
%              EF=-1: maximum number of iterations exceeded.
 
%  C) Quasi-Newton optimization algorithm using (BFGS)                  %
 
function [x,i,FunEval,EF]= Quasi_Newton (fun, x0, MaxIter, epsg, epsx) 
%   Variable Declaration 
 xi        = zeros(MaxIter+1,size(x0,2));
 xi(1,:)   = x0;
 Bi        = eye(size(x0,2));
 
%  CG algorithm
FunEval = 0;
EF = 0;
 
  for i = 1:MaxIter
 
      %Calculate Gradient around current point
      [GradOfU,Eval] =  Grad (fun, xi(i,:));
      FunEval        =  FunEval + Eval;       %Update function evaluation
 
      %Calculate search direction 
      di             = -Bi\GradOfU ;
 
      %Calculate Alfa via exact line search 
      [alfa, Eval]   =  LineSearchAlfa(fun,xi(i,:),di);      
      FunEval        =  FunEval + Eval;       %Update function evaluation
 
      %Calculate Next solution of X    
      xi(i+1,:)      =  xi(i,:) + (alfa*di)';
 
      % Calculate Gradient of X on i+1
      [Grad_Next, Eval] =  Grad (fun, xi(i+1,:));
      FunEval           =  FunEval + Eval;       %Update function evaluation
 
      %Calculate new Beta value using BFGS algorithm            
      Bi                =  BFGS(fun,GradOfU,Grad_Next,xi(i+1,:),xi(i,:), Bi);         
 
      % Calculate maximum acceptable Euclidean norm of the gradient
      if norm(Grad_Next,2) < epsg
          EF        = 1;
          break
      end
 
      % Calculate minimum relative change in the optimization variables
      E            =   xi(i+1,:)- xi(i,:);
      if norm(E,2) < epsx
          EF       = 2;
          break
      end
  end
  % report optimum solution
   x    = xi(i+1,:);
 
  if i == MaxIter
  % report Exit flag that MaxNum of iterations reach      
     EF =  -1;
  end
 
  % report MaxNum of iterations reach  
  Iter  = i;
 
end
 
%***********************************************************************%
% Broyden, Fletcher, Goldfarb and Shanno (BFGS) formula
%***********************************************************************%
function  B  = BFGS(fun,GradOfU,Grad_Next,Xi_next,Xi,Bi)
 
 % Calculate Si term
  si               =  Xi_next   - Xi;
 
 % Calculate Yi term
  yi               =  Grad_Next - GradOfU;
 %
 % BFGS formula (Broyden, Fletcher, Goldfarb and Shanno)
 %
  B   =  Bi - ((Bi*si'*si*Bi)/(si*Bi*si')) + ((yi*yi')/(yi'*si'));
 
end

References

Bonnans, J. F., Gilbert, J.Ch., Lemaréchal, C. and Sagastizábal, C.A. (2006), Numerical optimization, theoretical and numerical aspects. Second edition. Springer. ISBN 978-3-540-35445-1.
William C. Davidon, Variable Metric Method for Minimization, SIOPT Volume 1 Issue 1, Pages 1–17, 1991.
Fletcher, Roger (1987), Practical methods of optimization (2nd ed.), New York: John Wiley & Sons, ISBN 978-0-471-91547-8 .
Nocedal, Jorge & Wright, Stephen J. (1999). Numerical Optimization. Springer-Verlag. ISBN 0-387-98793-2.
Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007). "Section 10.9. Quasi-Newton or Variable Metric Methods in Multidimensions". Numerical Recipes: The Art of Scientific Computing (3rd ed.). New York: Cambridge University Press. ISBN 978-0-521-88068-8. http://apps.nrbook.com/empanel/index.html#pg=521.

Optimization: Algorithms, methods, and heuristics

Unconstrained nonlinear: Methods calling ...

... functions

Golden section search · Interpolation methods · Line search · Successive parabolic interpolation

... and gradients

Convergence	Trust region · Wolfe conditions

Quasi–Newton	BFGS and L-BFGS · DFP · Symmetric rank-one (SR1)

Other methods	Gauss–Newton · Gradient · Levenberg–Marquardt · Conjugate gradient

... and Hessians

Newton's method

Constrained nonlinear

General	Barrier methods · Penalty methods

Differentiable	Augmented Lagrangian methods · Sequential quadratic programming · Successive linear programming

Convex minimization

General

Interior point method · Reduced gradient (Frank–Wolfe) · Subgradient method · Cutting-plane method

Linear and
quadratic

Interior point	Ellipsoid method of Khachiyan · Projective algorithm of Karmarkar · Semidefinite programming

Basis-exchange	Simplex algorithm of Dantzig · Criss-cross algorithm · Principal pivoting algorithm of Lemke

Combinatorial

Paradigms

Approximation algorithm · Dynamic programming · Greedy algorithm · Integer programming (Branch & bound or cut)

Graph algorithms

Minimum spanning tree	Bellman–Ford algorithm · Borůvka · Dijkstra · Floyd–Warshall · Johnson · Kruskal

Network flows

Dinic · Edmonds–Karp · Ford–Fulkerson · Push-relabel maximum flow

Metaheuristics

Evolutionary algorithm · Hill climbing · Local search · Simulated annealing · Tabu search

Categories (Algorithms · Methods · Heuristics) · Software

Quasi-Newton method

Contents

Description of the method

Example algorithm using matlab

See also

References